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A new method to obtain a local parameterization for the exchange term in the many- 
\/~} • body electronic problem is presented. The approach amounts to the introduction 

of a coordinate dependent electron effective mass. Numerical results for metallic 
clusters in the jellium model are compared with other standard methods. 
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I. INTRODUCTION 



One of the major challenges in the density functional theory [|T]-§] is to improve the 
treatment of the exchange and correlation energy terms in inhomogeneous systems. Among 
the most vigorously pursued schemes is the inclusion of generalized gradient corrections to 
the local density approximation (LDA) One of the persistent problems related to the 

current formulations of LDA is the fact that the single-particle energies obtained in these 
methods have no direct relationship with the actual single-particle spectrum of the systems 
under consideration (with the exception of the ionization energies though). In particular, the 
gaps in semiconductors are severely underestimated || . It has been argued 0] that the origin 
of this discrepancy is the very nature of the local approximation to the density functional 
theory. This is one reason why, following an earlier argument ||, the electronic properties 
of several compounds were recently computed in an exact treatment of the exchange energy 
0. An improvement to the gap problem in semiconductors seems to be provided by the so 
called self-interaction method JT0[ . 
^ ■ Here we have chosen to explore another line of inquiry, inspired by the so called opti- 

mized effective potential (OEP) treatment of the exchange energy. This method was first 
introduced in atomic physics |TT| , p!2| and has been increasingly revisited lately ||13| . This 



recent activity generated some significant developments, for instance, a time dependent ex- 
tension of the OEP approach |14] . So far, most of the applications of the OEP have been 



devoted to the exchange only functionals. 

As in the OEP method, we shall consider the total energy of a many electron system 
E tot in the Hartree-Fock approximation only 
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where mo is the electron mass, e is its charge, Ri gives the position of nuclei with charge 
Zie, ^(r, er) the single-particle electron wave functions and a the spin variables. We shall 
use throughout this work the index h for labelling the hole (occupied) states and p for the 
particle (unoccupied) states. The indices h and p will stand for the corresponding quantum 
numbers of the single-particle states. The single-particle density is given by 

p a (r, r') = £ 9 h {r, Wh^, *), p{r, r>) = £ p a (r, r') . (2) 

h tr 

For the sake of simplicity we shall suppress the spin variables in the following formulas. 

The standard approach to self-consistently minimize E tot is to solve the Hartree-Fock 
equations. As customary, those are obtained by varying E tot with respect to ^* h {r) keeping 
the single-particle wave functions normalized, leading to 

flHF*fc(r) = (~ V 2 + V diT (r)) * h (r) ~ e 2 £ / dr' ^ r '^ h '^ * h ( r >) = e h * h (r) . 

\ ^TTIq ) h' \V V \ 

(3) 

The local (direct) part of the potential is given by 

Vdir(r) = V ions (r) + e 2 J dr' r^K • (4) 

In the OEP approximation the single-particle wave functions are the solutions of a local 
Schrodinger equation 

h 2 

- —V 2 ^ h {r) + K,Ep(r)*&(r) = eA(r) (5) 

and the potential Vqep(t) is determined so as to minimize the total energy of the system 

mm{E tot ({VoEr(r)})} =► J*™ = ■ (6) 

At this point a comment about the well known Slater approximation for the exchange energy 
is in order. We shall call Slater the approximation in which the single-particle wave functions 
obtained in the LDA with exchange only (LDAX) are used to compute E tot according to 
Eq. (|Tj) . Thus, it should be fairly obvious that total energy estimates in these three methods 
satisfy the following relation 
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Efl < £° EP < « ter • (7) 

In the same spirit, notice that one cannot vouch for the value of E tot computed in LDAX 
(i.e. consistently using the Slater prescription for the exchange energy to calculate E tot ) to 
be either an upper or lower bound estimate for the total energy. 

This paper is organized as follows: In the next section we present a generalized local 
approximation (GLA) and explain the naturalness of introducing a coordinate dependent 
effective mass. In Section III we describe two strategies for implementing the GLA and 
discuss in detail the case of spherically symmetric systems. In Section IV we present results 
for metallic clusters in the jellium approximation. There we compare different approximation 
schemes. We comment on some further possible extensions of GLA in Section V and present 
our conclusions in Section VI. 



II. THE GENERALIZED LOCAL APPROXIMATION 

In this section we shall present an extension of the OEP method, which we shall refer 
to as the generalized local approximation (GLA). This method is partially inspired by early 
attempts to treat the exchange term in a systematic way, put forward mostly in nuclear 
physics [p!5| -[T8[1. The basic idea of the GLA is to replace the OEP Schrodinger equation 



(^) by a generalized local Schrodinger equation with a coordinate dependent effective mass 
m eff (r) = m n(r) 

h 2 1 

- - V* A (r) + V{r)V h {r) = e h V h {r) . (8) 



2m n{r) 

The GLA local potential V(r) and the effective mass n(r) are now determined by following 
set of equations 

rmn{E tot ({V(r),n(r)})} =^ ^ = « ^ = 0. (9) 

To motivate that a coordinate dependent effective mass is a natural ansatz, we shall invoke 
two different arguments. The first one, formal and general in nature, is essentially a summary 



of a more comprehensive reasoning presented in Ref. [Tj|. The second argument is more 
physical, but will be presented by putting in perspective the different approximation schemes 
for the particular case of a jellium model for metallic clusters. 

Let us start by writing the Hartree-Fock equations (|[) in the form 

- — - V 2 ^(r) + / dr' U(r, r') % h {r') = e h ^ h (r) , (10) 

where the kernel U(r, r') contains the exchange potential which we are interested in, plus a 
direct term Vdh{r)5(r — r') as given by Eq. (f|). For convenience, we define the new space 
coordinates 

x — -(r + r') and s = r' — r (11) 
2 
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and change accordingly the kernel to 



U{x, s) = U{r, r') . (12) 

An approximate way to obtain a local equivalent of (|T0D is to expand the wave functions and 
the diagonal part of the kernel U in a Taylor series in s, retaining terms up to the second 
order in s 

Mr') - **(r) + V r * fc (r) ■ s + \ E (13) 
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and the kernel itself 

U(r + Is, s) » U(r, s) + 1 V r U(r, *) ■ * + \ £ > ( 14 ) 

where z and j label the axis and rj is the component of the vector r along the axis i. Direct 
insertion of the approximate Nt^r') and U(r + |a, s) into Eq. (|10D gives an equation of 
the very same structure as Eq. (||). Furthermore, we can explicitly write for the optimized 
potential 

V(r) = / ds U(r, s)+ 1 -Y,fds (15) 



and for the effective mass 



jj,(r) h 



^E dsU{r,s)siSj . (16) 



This set of equations is a good local approximation to the exchange term provided that 
the range of non-locality so in U(r, s) is smaller than the local wave length, or in general, 
smaller than the typical bulk characteristic length of the direct part, i.e. sq -C 27t//cf, where 
hp is the Fermi wave number. This approximation scheme was put forward by Frahn and 



Lemmer |[5] . A comprehensive overview on this and similar approaches is presented in Ref . 
0- 

One might be tempted to use Eq.flToD and flTED to construct V(r) and The problem 

is that often we cannot assume kpSQ to be small, particularly at surfaces. Unfortunately, 
a well defined systematic approximate scheme to construct V(r) and fi(r) is only known 



for non-local potentials in one-dimension [18|. Even for the simple spherically symmetric 



three dimensional case a satisfactory solution can be found only in some particular cases. 



In spite of this, the discussion presented above and the insight provided by Refs. pr^ , [16|JT8 
motivated us to consider the variational ansatz in Eq. @ as a very natural one. 

Another formal reason to introduce an effective electron mass in addition to an optimized 



local potential was pointed out by Austern | 17fl . The local and non-local wave functions 



differ in one rather subtle aspect. A node of a local wave function coincides always with 
an inflection point of that wave function, as one can easily see from the radial Schrodinger 
equation, if <f>(r) = then 0"(r) = as well. This is no longer true for non-local wave 
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functions and the only simple way to remedy this functional difference between local and 



non-local wave functions is the introduction of a position dependent effective mass ||16|| . The 
exact Hartree-Fock single particle wave functions, being solutions of a non-local equation, 
do not have the property that nodes correspond to inflection points. This is one reason why 
the trivial exact local effective potentials corresponding to wave functions with nodes have 
usually pole singularities. The remnants of these poles can be seen, for example, in the local 
effective potentials determined by Talman and Shadwick 
Based on these considerations one should expect that 

Eg < Ef^ < < £gf« , (17) 

since the GLA variational ansatz is more flexible and more suited for the exchange term 
than OEP. 

In the next sections we shall exemplify how the GLA method can be implemented on 
alkali atomic clusters |Tj|. Presently we shall treat the ionic background in the well known 
jellium approximation and consider only spherical closed shell clusters, when only the valence 
electrons are explicitly taken into account. One reason for choosing this system is that 
the electronic density of an alkali cluster is almost constant in the interior of the cluster. 
Inhomogeneities of the electronic distribution arise only because of the presence of a surface. 
The surface induces the natural falloff of the electronic density outside of the cluster, as well 



as Friedel-like oscillations ||19|| . The simplicity of the jellium approximation and the spherical 
symmetry of the closed shell clusters allows us to more easily single out the effects originating 
from inhomogeneities and the role of an effective mass. 

One of the appealing features of the OEP method for finite Coulomb systems is the 
correct asymptotic behaviour of the local potential, namely that for r — > oo, V(r) — > — e 2 jr. 
This feature is also characteristic of the GLA approach. Besides providing a better estimate 
for the total energy, which in itself might not be really a significant gain, the GLA method 
has an additional desirable feature: It generates a better approximation to the single-particle 
spectrum than the LDA and the OEP methods. 

Let us consider, for example, the Nag 2 cluster in various approximations. For more details 
see Ref. and Section IV. In the LDAX the width of the occupied band is A^ldax = 2.55 
eV, while in the Hartree-Fock approximation A^hf = 5.27 eV. The Fermi gas estimate gives 
Aepg = £>F/2mo, where pp is the Fermi momentum. For a Wigner Seitz radius r s = 4 a.u., 
which is approximately the value for bulk sodium, Aepg = 3.13 eV. The introduction of the 
self-interaction correction (SIC) in the exchange energy increases LDA band width from 
2.55 eV to 2.94 eV, which is closer to the Fermi gas estimate. This is an indication that 
the observed smaller occupied band width in clusters than in the bulk is due to electronic 
spill-out effects. The rather big discrepancy between the Hartree-Fock value and the Fermi 
gas estimate can be naturally attributed to an electron effective mass arising from the non- 



local Fock potential |I5]| . A naive estimate gives m e s/m = A£ Fg /Ae H F ~ 0.6. This value is 
actually very close to what we determine for this cluster. The reason for this is that in the 
cluster interior the electronic density is to a fair approximation equal to the jellium model 
bulk constant density. The occurrence of an electronic effective mass, different form the bare 



one, can be also interpreted as an energy dependence of the electron self-energy [21 . 

The trivial exact local potentials, for nodeless wave functions, determined from the exact 
solution of the Hartree-Fock equations for Na 92 system show that V-^ F ~ —9.12 eV and 



5 



~ —6.40 eV |2(J. These values can be interpreted either as an energy or as an angular 
momentum dependence of the trivial exact local effective potential. An angular momentum 
dependence of the effective local potential is in principle present and shall be discussed below. 
We claim, however, that the potential depth difference between V^ F and V^ F can be mainly 
accounted for by the energy dependence of the local effective potential. In other words, the 
effect can be interpreted by an electron effective mass smaller than the bare electron mass. 
In this model, for Nag2 the last occupied level lh {e\h — —3.38 eV) is very close in energy to 
the 3s level (e 3s = —3.50 eV), while the lowest occupied level Is has the energy e ls = —8.65 
eV. Had one attributed the energy difference Sxh — S\ s = 5.27 eV to an angular momentum 
dependence of the local effective potential, one would have great difficulty in explaining why 
the energy difference £3^ — E\ s = 5.15 eV is almost as big and not much closer to the Fermi 
gas estimate £F g = 3.13 eV. 

Before concluding this section, it is worthwhile to mention that in nuclear physics another 
approximation beyond the Slater prescription for the exchange energy has been suggested. 
The method, called density matrix expansion (DME), proved to be very successful for short 
range nuclear forces between fermions [ 22"f . In spirit, the DME method is a generalization 
of the traditional Slater approximation to inhomogeneous systems. The resulting local self- 
consistent equations for the single-particle wave functions are similar in structure to Eq.(^). 
The difference between DME and GLA is that in DME fi(r) and V(r) are self-consistently 
obtained from the eigenstates. Unfortunately, when applied to Coulomb systems the ensuing 
DME equations have inherent instabilities . 



III. NUMERICAL IMPLEMENTATION 

The basic novel idea of the present work is entirely contained in Eq.(|8]), as discussed in 
the previous section. Nonetheless, in order to have a useful generalized local approximation 
one needs an efficient minimization algorithm for the total energy E tot . This section is 
devoted to the discussion of two different strategies conceived to determine the best local 
potential and effective mass in the GLA approximation. 



A. Explicit parameterization 

In this first approach we shall represent the potential V(r) and effective mass /j,(r) as 
functions determined by a set of parameters {a^}. Hence, the problem of finding the best 
local potential is equivalent to finding the minimum of the functional E to t{{dk})- It is 
obvious that the quality of the GLA depends on the choice of the variational ansatz. In this 
section we shall first present the method without specifying any parameterization. We then 
proceed showing how this method is implemented for a spherically symmetric system. The 
numerical results shall be discussed in section IV. 

Recalling Eq.(||) the partial derivatives of the total energy with respect to {a^} can be 
written as 

9E M xfdrtZ&HnW+c*., (18) 



da k V J da 
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where d^h(r)/dak is the solution of the equation 

fir \d$h(r) (dH GLA de h \ 

where Hqla is the GLA Schrodinger operator. We stress that the conditions 

V(r) — > and ji(r) — > 1 for r — > oo (20) 

are satisfied throughout. 

Since the systems we consider in this paper are spherically symmetric one can write V(r) 
and fi(r) as functions of the radial coordinate r only. Thus, a suitable variational ansatz 
reads 

N 

V(r) = £ V k exp 



-N 



(r - r k )' 
2a? 



N 



and ju(r) = £ yUfcexp 



k=-N 



(r - r k f 



2a 2 



(21) 



where = afc. The sum over fc encompasses negative values in order to guarantee the 
correct behaviour near the origin, namely V'(0) = and //(0) = [p^| , with 14 = and 
jtife = Jl-k- Working with this ansatz implies that {a k } = ({Vk}, {Uk})- 

The method is implemented as follows: For a given starting V(r) and /x(r) in the form 
({Vk},{fik}) we solve Eq.(|). New amplitudes {V^} and {fik} are obtained by using the 
gradients given by Eq. ( jl8D in a suitable minimization algorithm. In particular, we used 
the molecular dynamics method described in Ref. |25|] and a simplified simulated annealing 
procedure in order to reach the condition dE tot jda k ~ and thus minimize E tot as given by 
Eq.©. 

B. Unconstrained minimization 

In this second approach we have applied the steepest descent method to find directly in 
coordinate representation the optimized local potential and the effective mass. At a first 
glance this approach is more attractive than the previous one, since it does not rely on a good 
variational ansatz. On the other hand, for practical use, an explicit parameterization allows 
one to reduce the space of minimization variables and obtain very efficiently reasonable 
solutions. However, if the system in question does exhibit a specific symmetry this method 
is particularly easy to implement. In the remaining of this section we discuss the particular 
case of spherical symmetric systems. We believe that this approach can also be implemented 
for other situations. 

In order to solve Eq. (|8|) numerically, it is convenient to represent the single-particle wave 
functions as 



* 4 (r) = Mr) Y hmi (r) = ^W) Y hmi (f) (22) 

where the index i labels states throughout the spectrum of the Hqla- After some straight- 
forward manipulations, Eq.(pl) can be written as 
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- (r) + (u(r) + l Ak±H _ Kr) e \ 0.( r) = o 



(23) 



where the energy was rescaled as 6j = 2m /h 2 Si and the potential Z7(r) is given by 



C/(r) 



2m 
IT' 



{5 









1 n"(r) n'(r) 

2 /x(r) r/x(r) 



(24) 



The normalization and completeness relations for the single-particle wave functions 0i( r ) 
read 



(25) 



£ 0fc(r')^(r)0 fc (r) = 5(r - r') . 



where J2k includes an integration over the continuous (unbound) spectrum of H GL ^. The 
functional variation of the total energy (|T|) can be brought to the following form 

5E tot = Jdr j<5f/(r) £(2fo + l)0 ft (r) Xh (r) - Spi(r) Y&h + l)e^(r) Xfe (r)| + c.c. . (26) 

The auxiliary functions Xh{ r ) are solutions of the following set of equations 



- xUr) 



jr, s . h{h + 1) / s 
U[r) H ^ M r K 



(27) 



where p is the single-particle density matrix as in Eq.(^), and huF refers to the radial 
part of the Hartree-Fock hamiltonian, Huf- On the right hand side of Eq . (|27D , [1 — p] 
is the projection operator outside of the occupied single-particle space and [1 — p]i h its 
radial component corresponding to the angular momentum 1^. In Eq. (p6|) the summation 
is over the occupied (hole) states h only The auxiliary single-particle functions Xhij') are 
orthogonal to the occupied single-particle states 



(0fcNxz) = / dr<f)* k (r)ii(r)xi{r) = 5 M 



(28) 



for all occupied states k and I. It is well known that an equivalent formulation of the 
Hartree-Fock equations is 



[1 - p}H nF p = , or [H UF } ph = [H nF } hp = or (1 - p)H UF ^ h (r) = 



(29) 



i.e. that all the particle-hole matrix elements of the Hartree-Fock single-particle Hamilto- 
nian vanish. Consequently, if the single-particle wave functions <fih( r ) are the solutions of 
the Hartree-Fock equations (||), then Xhij) = 0. One should not construe this as a statement 
that the Hartree-Fock minimum can be reached exactly in this way though, barring a pure 
coincidence. According to Eqs. @ and @ the single-particle wave functions have to be 
solutions of both GLA and HF equations at the same time, which in general is very unlikely. 
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The best one can hope for is that ( p9|) will be rather well satisfied for some U(r) and n(r). 
Since the total energy is bounded from below the existence of an optimized local potential 
and an effective mass is certain. However, the uniqueness of a minimum, i.e., the existence 
of only one global minimum, is not guaranteed. In principle one has the same problem in 
the exact HF case also. The vanishing of 5E to t upon variation of the local potential U(r) 
and effective mass /x(r) leads to the following equations 

J2(2h. + l)Mr)Xh(r) = , (30) 
h 

*E(2l h + l)e h <t> h {r) Xh {r) = . (31) 
h 

In the case when /i(r) = 1, Eq. (30) can be readily rewritten in the form of the equation for 
OEP derived by Talman and Shadwick |12|]. To implement the steepest descent method one 



has to change the local optimized potential and the effective mass according to the simple 
rules 

U(r) -> U(r) + SU(r) and fx(r) -> n(r) + 6fi(r) (32) 

where 

5U(r) = -2\ y £(2l h +l)<i> h (r)x h {r) 

h 

5fi{r) = 2A J2Wh + l)e h Mr)Xh(r) . (33) 
h 

The step size A has to be gauged with a certain care, so that Eq. leads to new single- 
particle wave functions corresponding to a lower total energy E tot . It is worth mentioning 
that the corrections to U(r) and /i(r) introduced above satisfy the following constraints 

/ dr fi(r) 6U(r) — and / dr fj,(r) 6fi(r) — , (34) 

which in particular imply that one cannot change the real effective local potential V(r) by 
a constant only, i.e. V(r) — > V(r) + const is not a possible change within this scheme. 

IV. RESULTS 

In our numerical study we used the jellium model to illustrate and compare the previously 
discussed approximation schemes. The model is defined by replacing V ions (r) as appears 
in Eq.(H), by the potential given by an uniform positive background charge density. We 
analyzed the following alkali clusters: Na 40 , Na 92 , Nai 38 and Nai 96 . For those, corresponding 
to electronic magic numbers, the spherical jellium approximation provides remarkably good 
results for the ground state properties and optical response The spherical geometry 



and the bulk Wigner-Seitz radius for sodium, r s = 4 a.u., determine the model entirely. In 
this section we examine the results of different schemes to solve the problem posed by Eq. ([!]) 
in the jellium approximation, namely Hartree-Fock, OEP, GLA and Slater. 
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Since the discussed schemes are approximations to the exchange potential, we first need 
to solve the jellium Hartree-Fock problem to have a standard to compare with. For this 
purpose we wrote a code that uses the same kind of algorithm as Ref. [2(1. Apparently the 
method used by Hansen and Nishioka |2(J is more accurate, but since we had no problems 
with obtaining converged results for cluster sizes up to Nai96, we did not improve further our 
code. As for the Slater functional, we used a method very similar to the one described in Ref. 
|27J , but with a series of refinements to ensure a higher numerical accuracy and increase the 
speed of computation. To obtain the solution of Eq. (|8|) and (§) we typically used the explicit 
parameterization method with iV = 20, • • ■ , 30, depending on cluster size and a = 1 A. By 
starting with a reasonable guess for /i(r) and V(r) we were always able to find a value for E tot 
very close to the Hartree-Fock value (see Tables I-III). Using the unconstrained minimization 
we could only marginally improve the minima obtained by explicit parameterization. In 
order to get a feeling of how well the GLA method works for different electronic densities we 
have varied the jellium density by a factor of two in both directions, i.e. smaller and higher 
densities (3 a.u. < r s < 5 a.u.). Again the results were always in very good agreement with 
the Hartree-Fock ones. 

We should mention, however, that for alkali clusters in the jellium model the GLA has 
likely several very close lying energy minima and often the corresponding U(r) and fi(r) differ 
quite considerably from each other. Furthermore, even for the states we have identified, we 
can only claim that our numerical solutions are very close to the actual minimum. The 
direct minimization procedures we have used, meet with considerable numerical accuracy 
problems close to a minimum and one can hardly improve on the quality of an approximate 
solution. The total energy of alkali clusters, as of any many-electron system as well, comes 
as a result of a strong cancellations of different contributions. Thus, a numerically accurate 
solution (relative accuracy better than pa 10~ 5 ) is very difficult to obtain. 

In Tables I-III we show the kinetic, Hartree, exchange, electron-ion and total energies 
for various sodium clusters in Slater, OEP, GLA and Hartree-Fock approximations. As one 
can see from these tables the GLA results are very close to the Hartree-Fock minimum and 
are of better quality than the OEP and Slater results. In Fig. 1 we compare the electron 
density profiles obtained in OEP, GLA and Hartree-Fock methods. Again, even though our 
solution might not exactly correspond to a local minimum, it is much closer to Hartree-Fock 
than the OEP solution. The local effective potential V(r) for the GLA and OEP cases are 
displayed in Fig. 2 and in Fig. 3 we show the renormalized effective potential U(r), see Eq. 
(23), and the effective mass /i(r). We have not imposed the "correct" asymptotic behavior 
-e 2 /r on the local potential as was done, for example, by Talman and Shadwick |12[]. Since 



at radii a few A larger than the jellium radius the electronic density is very small, any 
change of the potential or of the effective mass in this region has little effect on the hole 
single-particle wave functions. This is the origin of possible unexpected features in the local 
effective potential, such as non-monotonic behavior beyond the jellium radius. If unphysical 
characteristics appear in the local potential for large r at the early minimization stages, 
they are very difficult to correct. Such unwanted features can only significantly alter the 
particle (unoccupied) states. Whenever we encountered such problems, after obtaining an 
approximate local minimum in the GLA, we proceeded as follows: Using the parameterized 
minimization method we constrained the amplitudes {flk} and {Vk\ for gaussians inside 
the jellium radius to be fixed, varying only the others. Although these corrections have a 
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negligible effect on E to t } the effective potential tends slowly to acquire the correct form for 
large r. 

The single-particle spectrum we have obtained for the occupied states is in very good 
agreement with the HF one. The corresponding single-particle spectra for the occupied 
states in either LDAX, Slater or OEP methods is much more compressed, as is shown in 
Fig. 4. These differences arise because the effective mass in GLA is smaller than mo inside 
the cluster (m e fr ~ 0.6m ), as we have alluded to in Section II. 

The unoccupied single-particle states show very interesting features. In the Hartree-Fock 
approximation the level density of the unoccupied states is too small and the gap at the 
Fermi level is too large. The reason is well known: the particle states in this approximation 
are computed in the field of a complete screened positive charge. On the other hand, in the 
OEP method the level density of the unoccupied states is large and the gap at the Fermi 
level is small. The main reason for it is that the VoEp{ r ) exhibits the correct asymptotic 
behaviour as r — > oo. This causes the major difference between the OEP and the LDAX 
single-particle spectra. In LDAX one observes an overall upward shift of the spectrum, 
mainly due to the too sharp fall-off of the potential. One can reconcile almost perfectly 
the HF and the GLA particle spectra by rescaling the Hartree-Fock potential by a factor 
(N e — 1) /N e , where N e is the total number of electrons. In this way we force by hand the HF 
potential to have the "correct" —e 2 /r asymptotic behaviour for particle states. Notice that 
the hole states automatically have the correct asymptotic behaviour built in. The pleasant 
feature of the GLA approach is that one obtains not only a correct band width for the 
occupied states, but also a correct asymptotic behaviour for the effective potential of the 
particle states. As a result the gap at the Fermi level is larger in the GLA method than in 
the OEP method and smaller than in the HF approximation. 



V. FURTHER GENERALIZATIONS 

The GLA we have presented here is not yet the most general one. As we have mentioned 
in the introduction, one can extend it by considering an arbitrary angular momentum depen- 
dence of either the effective local potential or of the effective mass. (In general, an arbitrary 
angular dependence of the effective local potential leads strictly speaking to a non-local 
potential.) For example, one can consider a GLA Schrodinger equation of the following type 



h 2 



2m 



V M (r)Vtf fc (r) + V(r)V h (r) = e h V h (r) ■ 



(35) 



where M (r) is a symmetric tensor of rank two. One can chose this tensor in an appropriate 
manner such that the Schrodinger equation still preserves the spherical symmetry. For 
example, in the case of spherical systems one can choose 



M (r) 



1 b(r) 



where % and j are axis labels; and rewrite Eq. 



(36) 



as 
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U(r) + h{h + 1) (l + fc(r))-/x(r)e, 



(r) = . (37) 



In this case one can interpret the appearance of a tensor inverse effective mass as an angular 
momentum dependence of the effective local potential U(r). This ambiguity occurs only 
in the case of spherical symmetry. A simple analysis of the HF single-particle spectra of 
the clusters discussed in the previous section suggests that such a correction exists. The 
fact that the splitting between consecutive s-levels is smaller than between p-levels, which 
in its turn is smaller than the splitting between cHevels and so forth, gives room for such 
speculation. On the other hand, the total energy and the electronic density distribution of 
the alkali clusters we have considered here is very little affected by this type of correction. 
We have investigated to some extent this possibility as well, but not very thoroughly. 

One can consider some further generalizations, e.g. to introduce a pseudovector compo- 

nent (or equivalently an antisymmetric component of M (r)) of the inverse effective mass 
and/or also the introduction of a vector effective potential as well and still have a local 
Schrodinger equation which is a partial differential equation of at most second order. For 
example, a natural candidate for a pseudovector is the spin density. One can expect that 
in spin unsaturated systems a term with this symmetry might arise not only in the local 
effective potential but also in the inverse effective mass as well. One has to keep in mind that 
this type of correction violates time-reversal symmetry (which is violated in these systems 
anyway) and consider it with care. Whether there will be a reason to introduce a vector 
effective potential as well is still unclear at the present moment. The presence of an effec- 
tive vector potential might signal also the existence of nonvanishing currents in the ground 
state, since in such a case the GLA Hamiltonian would not automatically be invariant under 
the transformation p — > —p, which will occur for example if the system is in an external 
magnetic field. 



VI. CONCLUSIONS 

Our analysis suggests a possible way to generalize the popular LDA treatment of Coulomb 
systems by enlarging the class of considered functionals. There is no a priori or fundamental 
reason why one should consider a total energy functional of the Kohn-Sham type only and 
not allow for terms, in which the inverse mass is replaced, e.g., by a density dependent 
function. One can find typical examples of such functionals in nuclear physics [f22j . By con- 
sidering such type of generalized local energy density functionals one can improve the quality 
of the description of inhomogeneous systems (total energy, electronic density distribution) 
and at the same time achieve a significant improvement of the single-particle spectrum as 
well, which is a long standing unsolved issue in LDA. We would like to remind the reader of 
one potential problem with using the present approach in infinite homogeneous systems. It 
is well known that the electron self-energy in the HF approximation for a pure Coulomb in- 
teraction between electrons has a logarithmic singularity at the Fermi level PTJ . In principle 



that should prevent us from introducing an effective mass the way we have described here. 
On the other hand the Coulomb interaction between electrons is strongly renormalized in 



medium, see e.g. Ref. [28], and this fact in particular leads to Coulomb screening. (In its 
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simplest classical form that is the Thomas-Fermi screening at large distances.) This should 
be sufficient to remove such singularities of the electron self-energy and thus lends support 
to a GLA approach. 
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TABLES 



TABLE I. Results for Nag2- Note that in LDAX and Slater the entries for the total kinetic 
energy £kim Hartree energy i^Hartreej total electron-ion interaction energy -Ecouiomb are identical 
(since these quantities have identical functional forms) and the only differences are in Fock or 
exchange energy -Epock and total energy of the cluster E tota \ respectively. Energies are in eV. 



LDAX Slater OEP GLA HF 

160.80 160.80 161.59 161.94 161.73 

^Hartree 7,560.71 7,560.71 7,565.32 7,565.09 7,564.02 

^Coulomb -15,210.17 -15,210.17 -15,214.67 -15,214.47 -15,213.37 

£ Fock -269.39 -277.28 -278.38 -279.18 -279.04 

£ total -188.41 -196.30 -196.50 -196.97 -197.01 



TABLE II. Results for Nai38- See caption for Table I for notation and additional remarks. 



LDAX Slater OEP GLA HF 

242.51 242.51 243.30 243.68 243.66 

^Hartree 14,889.18 14,889.18 14,902.39 14,902.44 14,901.75 

^Coulomb -28,926.61 -28,926.61 -29,939.80 -29,939.84 -29,939.15 

£ Fock -406.51 -415.98 -416.91 -417.85 -417.94 

£totai -268.37 -277.84 -277.96 -278.51 -278.61 



TABLE III. Results for Naig6- See caption for Table I for notation and additional remarks. 



LDAX Slater OEP GLA HF 

Fh~ 346.40 346.40 348.04 348.54 347.73 

^Hartree 26,776.91 26,776.91 26,791.67 26,792.13 26,781.13 

^Coulomb -53,764.64 -53,764.64 53,779.52 -53,779.76 -53,768.72 

E Fock -580.21 -591.13 -592.97 -594.33 -593.73 

Etatei -365.16 -376.08 -376.41 -377.05 -377.21 
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FIGURES 



FIG. 1. Electronic densities (in A~ 3 ) as a function or the radial coordinate r (in A) for 
Nag2 and Nai32 in the Hartree-Fock, GLA (solid line) , OEP (dashed line) and LDAX (dot-dashed 
line) approximations respectively. On this scale the GLA and Hartree-Fock densities are visually 
indistinguishable. 

FIG. 2. The optimized effective potential V(r) (in eV) for Nag2 in the GLA (solid line) and 
OEP (dashed line) approximations as a function of the radial coordinate r (in A). 

FIG. 3. For Nag2 the effective mass /x(r) (in dimensionless units) is displayed in the top panel. 
The renormalized optimized effective potential U(r) (in A~ 2 ), see Eq. (23), for the GLA (solid 
line) and OEP (dashed line) approximations as a function of the radial coordinate r (in A) is shown 
in the bottom panel. 

FIG. 4. Single-particle level scheme for Nag2 (in eV) in different approximation schemes. 
Solid lines correspond to the occupied spectrum and the dashed lines to several lowest unoccupied 
levels. For further details see text. 
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